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Abstract The effect on the convergence of warm start¬ 
ing the projected Gauss-Seidel solver for nonsmooth 
discrete element simulation of granular matter are in¬ 
vestigated. It is found that the computational perfor¬ 
mance can be increased by a factor 2 to 5. 
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1 Introduction 

In simulations of granular matter using the nonsmooth 
discrete element method (NDEM) [Il[2j[3] the computa¬ 
tional time is dominated by the solve stage, where the 
contact forces and velocity updates are computed. Con¬ 
ventionally this involves solving a mixed complementar¬ 
ity problem or a quasi-optimization problem that arises 
from implicit integration of the rigid multibody equa¬ 
tions of motion in conjunction with set-valued contact 
laws and impulse laws, usually the Signor ini-Coulomb 
law and Newton impulse law. The computational prop¬ 
erties of the solution algorithms for these problems are 
largely open questions, lacking general proof of exis¬ 
tence and uniqueness of solutions as well as of gen¬ 
eral proof of convergence and numerical stability [4]. 
The projected Gauss-Seidel (PCS) algorithm is widely 
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used. The popularity of PCS is likely due to having low 
computational cost per iteration, small memory foot¬ 
print and produce smooth distribution of errors that 
favour stable simulation. In many cases PGS require 
few iterations to identify the active set of constraints. 
This make PGS a natural choice for fast simulations of 
large-scale rigid multibody systems with frictional con¬ 
tacts. The asymptotic convergence, however, is slow. 
The PGS algorithm solves each local two-body contact 
problem accurately but approaches to the global solu¬ 
tion in a diffusive manner with iterations. This limit 
the practical use of PGS for simulations of high accu¬ 
racy. The residual error appear as artificial elasticity 
[5], with an effective sound velocity vpqs = 
where Na is the number of iterations, d is the parti¬ 
cle size and At is the timestep. Accurate resolution of 
the impulse propagation in stiff materials thus require 
large number of iterations or small timestep. The re¬ 
quired number of iterations for a given error tolerance 
increase with the size of the contact network, particu¬ 
larly with the number of contacts in direction of grav¬ 
ity or applied stress. It may, however, saturate by an 
arching phenomena analogous to the Janssen’s law for 
silos m- The PGS algorithm is parallelizable, despite 
many many authors claim of the opposite, for hardware 
with distributed memory using domain decomposition 
methods m\- 

Warm starting means to start the PGS algorithm 
with an initial guess, Aq", that presumably is closer to 
the exact solution, A , than starting with the nominal 
choice of Ao = 0. The idea, illustrated in Fig. [U is that 
the warm started PGS reach an approximate solution, 
A/c/, with fewer iterations than the solution, A/^, start¬ 
ing from nominal value. In other words, |A — A^/| < 
|A — A/c| < £ with k' < k and error tolerance e. The 
effective increase in convergence should be most sig- 
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Fig. 1 Illustration of improved convergence by warm start¬ 
ing. 

nificant for static or nearly static configurations. For 
rapid granular flows the solution change rapidly with 
time and no or little effect on convergence is expected. 
There have been several reports on improved conver¬ 
gence by using warm starting plISlI^ fTOlfTTlfT^ but to 
the best of our knowledge no quantitative analysis has 
previously been presented. 

2 PGS for nonsmooth discrete element 
simulation 

The mixed complementarity problem (MCP) for com¬ 
puting the update of the velocity from Vom = v{t — At) 
to V = v{t) and the Lagrange multiplier A of the con¬ 
tact constraints and take the form 


'M -G"^' 


V 


P 

G S 


A 


.q. 


A(“) eC^(A(“)) , a = l,2,...,Ve (2) 

where M is the mass matrix and G the Jacobian of con¬ 
tact constraints. The contact force, G^A, is restricted 
by a friction cone condition that we represent G 
where a indexes the contacts. The diagonal 
perturbation U regularize the problem and allow mod¬ 
eling of contact elasticity. The vectors p and q on the 
right hand side depend on particle inertia, external force 
and constraint violations on position and velocity level. 
As friction cone condition we use the Signorini-Coulomb 
law including 0 < and | with 

the friction coefficient /ig for each contact a divided 
in one normal (n) and two tangential (t) components. 
The constraint forces act to prevent contact overlap, 
gf < 0, and contact sliding, GtV = 0. Similarly, rolling 
resistance (r) is imposed by a constraint G^v = 0 with 
a Coulomb like law: |Ar^^| < where 

r* is the effective radius. See Appendix A for further 
details. For a system with Np particles represented as 


rigid bodies and Nc contacts with normal and tangen¬ 
tial force and rolling and twisting resistance the vec¬ 
tors and matrices in Eq. m have the following dimen¬ 
sions dim(M) = 6Np x 6A^p, dim(G) = 6Nc x Np, 
dim(i;) = dim(p) = 6Np dim(A) = dim(gr) = QNc. The 
matrices are however very sparse. M and U are block 
diagonal and G is block sparse. The blocks have dimen¬ 
sion 6x6. The main steps of the PGS iteration are 

(3) 

4-projcyAi^“\) (4) 

Vfe+i = + (5) 

with iteration index k = 0,1,2,..., A^it — 1, change in 
multiplier AX^j^^ = — X^j^^ and residual 

^1°"^ = ^P(c)-q(c) = G(c,)V/c-q('^) 

( 6 ) 

where v/. = M“^p + and D is the block 

diagonal part of the Schur complement matrix S = 
GM“^G^ + U. The details of the vectors p and q 
depend on the stepping scheme and constraint stabi¬ 
lization method. When integrating with flx timesteps 
At using the SPOOK stepper m one has p = Mvom + 

Z\tfext: with smooth external forces fext, and q = (q^, qf, q^^ 
with q„ = -{4:/At)Yg + TGnVoid, qt = 0 and q^ = 0. 

The projection A^^^ t— proj^^ is made by simply 

clamping to the friction or rolling resistance limit 
if exceeded. After stepping the velocities and positions 
an impact stage follows. This include solving a MCP 
similar to Eq. o but with the Newton impact law, 
G^"^v+ = —replacing the normal constraints 
for the contacts with normal velocity larger than an im¬ 
pact velocity threshold The remaining constraints 
are maintained by imposing Gv+ = 0. An algorithm of 
NDEM simulation with PGS is given in Appendix A 
together with details on the Jacobians and relation be¬ 
tween the solver parameters and material parameters. 


By default the PGS algorithm is initialized with Xq^^ = 
0 . We refer to this as cold starting. In a stationary state 
the contact force G^A is constant in time. In a nearly 
stationary state we expect the multipliers to remain 
almost constant between two timestep. Therefore it is 
reasonable to use the solution from last timestep as 
an initial guess, A(t) ^ X{t — At). We use a fraction 
yd = 0.85 of the solution from last timestep 


Ao(t) = (3XN,,{t - At) 


3 PGS warm starting 


( 7 ) 
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It is important to also apply the corresponding impulse 
to the particles and update the velocity 

Vo = M“^p + M“^G'^Ao (8) 

such it become consistent with the initial guess for the 
multiplier. We refer to warm starting based on the last 
solution as history based warm starting. For any new 
contact we set = 0. Warm starting is not applied 
at the impact stage and we assume that the contact net¬ 
work is not fundamentally rearranged by the impacts 
and use the solution from last timestep despite the oc¬ 
currence of impacts. 

An alternative method for warm starting a nearly 
stationary state is to estimate each local contact force 
and assign this to the contact multipliers. When us¬ 
ing regularized NDEM the local contact force can be 
estimated from the overlaps and relative contact veloc¬ 
ities much as in conventional smooth DEM. We refer to 
this approach as model based warm starting. Eor nor¬ 
mal forces we use the Hertz contact law /n = k^gn^'^^ 
with overlap function and based on fn ~ G^X^/At 
we estimate 


Table 1 Main material and simulation parameters 


Notation 

Value 

Comment 

[d, d2] 

[13, 10] m 

bi-disperse particle diameter 

P 

3700 kg/m^ 

particle mass density 

E 

6 MPa^ 

Young’s modulus 

e 

0.18 " 

restitution coefficient 

Ms 

0.91 

surface friction coefficient 

Mr 

0.32 

rolling resistance coefficient 

At 

5 ms 

timestep 

'^imp 

0.05 m/s 

impact threshold 



iVit=10 iVit=50 Vi,= 500 Vir=10 Vir=5 


Fig. 2 Samples of five columns simulated, from left to right, 
with cold starting Nit = 10, 50 and 500, and history based 
warm starting = 10 and 5. 




5/4 

(a) 


(9) 


4.1 Column 


Similarly the regularized tangent friction force and rolling 
resistance force can be estimated via the Rayleigh dis¬ 
sipation functions to 


At,o 

Ar,0 


7t-iz\i(Gtv)TGt 

77Mt(GrV)'^Gr 


(a) 


Note that the friction and rolling resistance should, 
if large, be clamped to obey the conditions |A 
and 


Particles of diameter d = 13 mm are initiated on top of 
another with zero overlap. The system compress slightly 
under its weight. The simulation is run until the ID 
column have come to rest. Sample images from simula¬ 
tion with and without warm starting and for different 
number of iterations are shown in Fig. [2j Warm start¬ 
ing clearly improve the convergence. To make a quan¬ 
titative convergence analysis we study the deviation of 
t I — the simulated column height, from the theoretical 
height, /, computed using the Hertz contact law 


( 10 ) 

( 11 ) 


= 


I - ^Nit 


4 Numerical experiments 

Numerical simulation of different systems were performed 
to analyse the effect of warm starting on the conver¬ 
gence of NDEM simulations. The following systems were 
studied: a ID column, formation of a pile, dense flow in 
a rotating drum and a triaxial shear cell. The main ma¬ 
terial and simulation parameters are listed in Table [H 
The method was implemented in the software AgX Dy¬ 
namics m in the module for NDEM simulation with 
optimized data structures and support for collision de¬ 
tection and PGS using parallel computing on multicore 
processors. The simulations were run on a desktop com¬ 
puter with Intel(R) Core(TM) i7 CPU, 2.8 GHz, 8 GB 
RAM on a Windows 64 bit system. Videos from simula- 


l 


( 12 ) 


A series of simulations are run with number of par¬ 
ticles, Vp, ranging from 5 and 100, number of itera¬ 
tions, Vit, ranging from 10 to 500. The required number 
of iterations, to reach a solution with error toler¬ 
ance s\ = 0.1%, 1% and 5% are presented in Eig. [3l 
It scales almost linearly with the number of particles 
and increase with decreasing error tolerance s\. History 
based warm starting is on average three times as effi¬ 
cient as cold starting. Also model based warm starting 
improve the convergence at low error tolerance. The 
performance gain from model based warm starting de¬ 
crease with increasing error tolerance and for s\ = 5% 
model based warm starting require twice as many iter¬ 
ations as cold starting. Eigure 0 ] show the evolution of 


tions are found at http: //umit. cs . umu. se/granular/wariiteaii^baimgei'sidual, see Eq. for the normal force con- 
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Fig. 3 The required number of iterations for a ID column 
simulation for error tolerance eir> = 0.1% (top), 1% (middle) 
and 5% (bottom) depending on the number of particles and 
warm starting method. 


straint during a PGS solve for a column with Np = 25. 
The convergence rates are similar but warm starting 
clearly has the advantage of starting closer to the solu¬ 
tion. 



Fig. 4 The evolution of the mean normal force residual dur¬ 
ing a PGS solve for a Np = 25 column using cold starting 
and history based warm starting. 


4.2 Pile formation 

A pile is formed by continuously emitting particles of 
diameter d = 13 mm from a 3d wide source placed 20d 
above a ground plane. The number of particles in the 
pile is Np = 3363. Again we use the relative height, e\ 
in Eq. m, as error measurement. The reference height 
about 15d is measured from a pile constructed using 
small time-step At = 0.2 ms and A^it = 500. Pile for¬ 
mation is then simulated using time-step At = 5 ms 
for different number of iterations and warm starting 
methods. Sample images from simulations are shown 
in Fig. [5l The angle of repose is an alternative mea- 



Fig. 5 Samples from simulations of pile formation. From left 
to right is the cold started pile {At = 5 ms, Nt = 50), a 
reference pile {At = 0.2 ms. Nit = 500) and a warm started 
pile {At = 5 ms, A^it = 50) 


sure but was found to give less precise result. The pile 
height is measured 10 s after the last emitted parti¬ 
cle has come to approximate rest. Simulations are run 
with warm starting applied both to normal forces, fric¬ 
tion and rolling resistance and to normal forces only. 
The historical warm starting is tested with and with¬ 
out the velocity update associated with the warm start 
in Eq. ©• The required number of iterations for a given 
error threshold are given in Fig. [6l With few iterations 
the piles experience artificial compression and contact 
sliding such that the pile gradually melt down to a singe 
particle layer. The pile stability increase with the num¬ 
ber of iterations. History based warm starting, applied 
to both normals, friction and rolling resistance, give 
the best result and require roughly half the number of 
iterations of cold starting. If the warm start velocity 
is not applied the result is worse than cold starting. 
Model based warm starting is only marginally better 
than cold starting and is from further experiments here 
on excluded. Applying warm starting to the normal 
constraints only does not improve the convergence sig¬ 
nificantly. The convergence is also analysed by studying 
the evolution of the Lagrange multiplier and the resid¬ 
ual. The relative error of the normal force multiplier is 
computed as 




\ i^S'i / 


(13) 


The evolution of during a solve of a stationary pile 
is shown in Fig. [71 The multiplier error for history based 
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for normal force constraints. This is consistent with the 
faster melting of the piles simulated with cold starting. 



Fig. 8 The mean residual dependency on the number of iter¬ 
ations when simulating a resting pile for 1 s using cold starting 
and history based warm starting. 


Fig. 6 The required number of iterations versus pile height 
error for different warm starting methods. 


warm starting is roughly five times smaller than for cold 
starting and remain more accurate indefinitely. A more 



Fig. 7 The evolution of the relative multiplier during a PGS 
solve for a resting pile using cold starting and history based 
warm starting. 


careful analysis can be made by studying the evolution 
of the residual, defined in Eq. and how it is dis¬ 
tributed over the constraints. To get comparable states 
a stationary pile is prepared by using 500 iterations 
from which the cold and warm started simulations are 
started and run for 1 s before the measurement. The 
evolution of the mean residual during a PGS solve is 
shown in Fig. [8l The convergence rates are similar but 
the initial lead of history based warm starting over cold 
starting by roughly a factor 5 remains throughout the 
500 PGS iterations. Gomparing the residual histograms 
from using cold and warm starting in Fig. [9] it is clear 
that the solutions differs primarily in the errors for the 
rolling resistance and friction constraints and less so 




residual 


Fig. 9 The residual distribution for a resting pile after 1 s 
using Nit = 100 iterations, cold starting (top) and history 
based warm starting (bottom). 


4.3 Rotating drum 

A cylindrical drum with diameter D = 40d and width 
w = 7d is rotated with angular velocity i? = 0.25 
rad/s. This corresponds to the Fronde number Fr = 
f2g ^ 10“^ which corresponds to the dense rolling 
flow regime. A nearly stationary flow of A^p = 4864 
particles with bi-disperse size distribution d and ^ 2 . At 
this low Fronde number a large plug-zone is developed 
where particles co-rotate rigidly with the drum. A con¬ 
vergence analysis is made of the plug-zone number frac¬ 
tion, A^piug/A^p, and the dynamic angle of repose, 0'. 
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These are measured for different number of iterations 
on a flow averaged over 2 s for cold starting and histor¬ 
ical warm starting. The sample trajectories in Fig. [10] 
illustrate the general trend that the dynamic angle of 
repose and the size of the plug zone decrease with de¬ 
creasing number of iterations but less so using warm 
starting. The normalized particle flow velocity relative 



Fig. 10 A sample of particle trajectories from simulation of a 
rotating drum with i? = 0.25 rad/s, At = 5 ms and Nit = 10 
(left), Nit = 500 (middle) and warm starting Nit = 10 (right). 


the plug flow is computed 'z;* = |v* — F x Q\/RQ and 
sample plots are shown in Fig. [TTl As threshold for the 
plug zone flow we set < 0.15, which is fulfilled by 
^phig/^p ~ ^ particles where the variations 

reflect the slightly pulsating nature of the flow, due to 
sequential onset of avalanches. The plug zone fraction 
number error is defined 


^plug — 


N500 _7vy 

^^plug plug 

Nr. 


(14) 


and the relation to the required number of iterations is 
found in Fig. [T2j The warm starting solution approach 
the solution faster but seems to have larger variations 
at high iteration numbers. The convergence analysis of 



Fig. 11 A sample of cross-section flow from a simulation of a 
rotating drum with i? = 0.25 rad/s, At = 5 ms and Ait = 10 
(left), Ait = 500 (middle) and warm starting Ait = 10 (right). 
The colour coding show the particle velocities relative to rigid 
co-motion with the drum. 



Fig. 12 The convergence of the plug zone fraction number 
for cold starting and history based warm starting. 



Fig. 13 The dynamic angle of repose as function of number 
of iterations for cold starting and history based warm start¬ 
ing. 


4.4 Triaxial shear 

The triaxial shear test is constructed by six dynamic 
rigid walls of mass 100 kg each that are driven with 
prismatic motors to apply a specific stress = fi/Ai, 
where Ai is the cross-section area and fi the applied 
motor force in the coordinate direction i = see 

Fig. HH First, a hydrostatic pressure of = 100 Pa is 



Fig. 14 Sample image from the triaxial test. 


the dynamic angle of repose also show that warm start¬ 
ing converges faster although to a slightly higher angle 
^w^piug ~ compared to = 48°, see Fig. [131 The 
dynamic angle of repose is measured as the displace¬ 
ment of the material centre of mass from the 2 ;-axis 
which is more robust than tracking the surface. 


applied on all sides. Then the top and bottom walls are 
driven inwards at 0.01 m/s by regulating and main¬ 
taining a constant side wall pressure at Cx = o-y = cr^. 
At some critical deviator stress cr° — the material 
fail to sustain further increase in stress and starts to 
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shear indefinitely. The transition is more or less sharp 
depending on the initial packing ratio, hydrostatic pres¬ 
sure and applied shear rate. In this test the Young’s 
modulus is set to the stiffer value of ^ = 60 MPa to get 
a sharper transition between compression and shear. 
The critical axial stress is computed as the averaged 
(Jz in the shear phase between lateral strain 5 = 10% 
to £ = 25%. The critical stress deviator depending on 
the number of iterations for cold starting and history 
based warm starting is shown in Fig. [151 Both curves 
converge to about 1 ± 0.2 kPa. With warm starting the 
stress levels out at Yit > 200 while cold starting re¬ 
quire Yit = 1000. Sample curves of the stress deviator 



Fig. 15 Critical yield stress as function of number of itera¬ 
tions for cold starting and history based warm starting. 

as function of lateral strain are shown in Fig. [161 These 
confirm the faster convergence when warm starting but 
also show higher stress fluctuations in the shear phase. 
Whether this is an artefact of the warm starting or an 
actual feature of the triaxial test has not been pursued. 



Fig. 16 Sample stress curves in triaxial test for 100 and 1000 
iterations. 


5 Application example 

The effect of using warm starting in practical simula¬ 
tion applications is illustrated with two examples. The 
first example is part of a balling drum circuit used in 
ore pelletizing systems m, see Fig.[T71 Simulations are 


used for the purpose of process control and for find¬ 
ing a design of the drum outlet that maximizes an 
evenly distributed throughput on a roller sieve where 
material is size distributed. Three distinctive subsys¬ 
tems with different dynamics can be identified. Firstly, 
there is the drum with an almost stationary flow. Sec¬ 
ondly, material is distributed onto quasistationary piles 
on a wide-belt conveyor. Thirdly, the particles disperse 
over a roller sieve with increasing gap size downwards 
to achieve a size separation. The design problem is fore 
mostly a geometric flow problem and the material dis¬ 
tribution need to be computed with sufficient accuracy. 
We assume 5 % is a required accuracy for dynamic and 
static angle of repose. From Fig. [13] we estimate that 
warm starting is roughly three times more computa¬ 
tionally efficient in computing the drum flow and, ac¬ 
cording to Fig. [6] twice as efficient for pile formation 
on the conveyor. The flow on the roller sieve is more 
disperse and collision dominated requiring only few it¬ 
erations (Yp < 25) and it can be expected that warm 
and cold starting are equally efficient. The overall com¬ 
putational speed-up by applying warm starting is thus 
estimated to a factor no larger than 2. 

The second example is an excavator. A rectangular 
trench is filled with roughly 10^ spherical particles of 
uniform size distribution between 25 and 100 mm and 
particle mass density of 2500 kg/m^. The excavator is 
modeled as a rigid multibody system of total mass 50 
ton divided in 10 bodies, 8 joints and 3 linear actuators 
(hydraulic cylinders) and one rotational motor. The full 
system of granulars and vehicle take the mathematical 
form of Eq. o and is solved using a split solver where 
the vehicle part is solved using a direct block-sparse piv¬ 
oting method [15] and the granular material with a PGS 
solver as described in this paper. Simulations were run 
with time-step h = 2.5 ms, which allow for a low num¬ 
ber of iterations. The machine perform an excavation 
cycle by a pre-programmed control signal to the actu¬ 
ators. The resulting actuator forces are measured and 
these include the back reaction from the resistance and 
inertia of the granular material. Two simulations, with 
and without warm starting, are run with Yit = 25. Sam¬ 
ple images from the simulations are shown in Fig. [18] 
Observe the difference in height surface of the granu¬ 
lar material due to artificial compression and frictional 
slippage due to numerical errors in the PGS solve. The 
undisturbed height in the two simulations differ by 10 
% and volume of displaced material differ by at least 
30 %. The difference in granular dynamics also affect 
the measured force response. The force trajectory of the 
middle actuator is provided in Fig. [19] In the phase be¬ 
tween 8 — 10 s, when the bucket is dragged through the 
material the force when using warm starting is almost 
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50% larger because more material is set in motion and 
stronger resistance to shear motion. Whether h = 2.5 
ms, A^it =25 and the improvement by using warm start¬ 
ing give sufficiently accurate force response depend on 
the intended use of the data and require further conver¬ 
gence analysis. On a desktop computeiH with the given 
NDEM settings the computational time is roughly 100 
s per realtime second. 



Fig. 17 A balling drum circuit with granular material in 
different states. 


6 Conclusions 

The convergence of the projected Gauss-Seidel algo¬ 
rithm for NDEM simulation is increased by warm start¬ 
ing with the solution from previous time-step. The com¬ 
putational speed-up by warm starting is demonstrated 
to be about 2 — 5 for a wide range of systems including 
pile formation, granular drum flow and triaxial shear. 
An examination of the residual distribution show that 
convergence improvement primarily improve on the ve¬ 
locity constraints - friction and rolling resistance - and 
less so on the normal force constraints. Warm starting 
the Lagrange multiplier based on 85 % of the value from 
last time-step was found to give best results. Warm 
starting based on an explicit contact force model give 
only marginal speed-up, for example 20 % for a pile for¬ 
mation. This is not surprising since the damping coef- 
flcients in the dissipation models for sliding and rolling 
are not physics based and can only predict the value of 
contact forces in slide mode but not of stick mode inside 
the friction and rolling resistance limits. Eor materials 

^ Performance measurement are made on a desktop com¬ 
puter with Intel(R) Core(TM) Xeon X5690, 3.46 GHz, 48 GB 
RAM on a Linux 64 bit system. 



Fig. 18 An excavator digging in trench with 10^ particles, 
h = 2.5ms, Ait = 25 and using cold starting (top) and warm 
starting (bottom). The colour codes the particle height with 
red to blue ranging from 0 m to — 2 m. Gray particles are 
above 0 m. 



Fig. 19 The force trajectory of the middle link pistons of 
the excavator while digging with Ait = 25 using cold starting 
and warm starting. 


shearing under high stress, compared to the stress pro¬ 
duced by the materials own weight, warm starting show 
larger stress fluctuations than without warm starting. 
Whether this is an artefact or correct behaviour emerg¬ 
ing after further iterations has not been established. 
A more in depth analysis of systems under large stress 
should be made considering also alternative size of time 
step, shear rate, hydrostatic load stress and particle 
stiffness. 
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Appendix 


A. Simulation algorithm 


The algorithm for simulating a system of granular ma¬ 
terial using NDEM with PGS solver with warm starting 
is given in Algorithm [H Based on the Hertz contact law, 


Algorithm 1 NDEM simulation with warm started 
PGS solver 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 


> Time stepping 


> impacts 


set constants and parameters 
initial state: (xo,vo) 
for 2 = 0, 1, 2,, tjAt do 
contact detection 
compute g, G, D 
impact stage PGS solve (v^, A^) 

compute + TnGnV+ 

pre-step v = v+ + AtM“^fext 
A/cq = 0 or warm start A/c^ 
warm-step v = v + 

-> PGS solve for continuous contacts- 

for k — 1,. .. , Nit and while criteria{r) do 
for each contact a = 0, 1,.. ., Nc — 1 do 
for each constraint n of contact a do 

(«) _ _a(«) I 

- \(«) _ \(«) 

^n,k ^n,k ^n,k — l 

v = v + 

end for 

end for 

end for 

Vi+i = V > velocity update 

Xi+i = Xi Atvi-i-i > position update 

end for 


. («) _ . («) 
'^n,k — '^n,k- 


> residual 
> multiplier 
[> project 


each contact a between body a and b add contributions 
to the constraint vector and normal and friction Jaco- 


5(a) 

= nfa)(Xa4 

-d7- 

- Xf, - d 

b ) 


9(a) 

= <^7 ,eH = 5/4 



G(a) 

^na 

= ^^9(a) ' 

T 

-(di“^ X 

n(a 

g(«) 

= ^M(aC 

T 

(di“' X 

n(a))^ 

gW 

^ta 


r 

-(dS"> X t<"> 

)T- 




-(d‘“' X t<“> 

)7 


^tb 


NWT 

(dl”> X 




_ 2 

(d'"> X 





0lx3 

J«)T 

•'1 

0lx3 — 

j-WT] 

gW 

^ra 

= 

0lx3 

J«)T 

•'2 

0lx3 — 

^2 



_0ix3 n(«)T 

0ix3 

g(") 


0lx3 

^1 

0lx3 

, (a) 

T " 

= 

0lx3 

^2 

0lx3 

.(a) 

^2 

T 



0lx3 

_n(«)T 

0ix3 n(“)T 


(15) 


where and are the positions of the contact 
point a relative to the particle positions and x^. 
The orthonormal contact normal and tangent vectors 
are and 

The diagonal matrices and Schur complement ma¬ 
trix D are 

_ 4 . 

— 


At^ 1 + 4^ 

- -^^2NcX2Nc 

y — 

- -^^3NcX3Nc 


^NcXNc 


(16) 


rn = 


1 


-lATeXATe 

D = GM-^G^ + Z' 


1 m AZn ■ 

^ At 


The mapping between MGP parameters and material 
parameters are 

5n = e^/kn = 3eu{l - iy‘^)/EV^ 

Tn = max(nsZ\t,5n/7n) (17) 

7n ^ = knc/el 


where r* = (r^ +r 5 )/rar 5 is the effective radius and we 
use 7t = 7r = 10“ 


- - - ris = 2. 
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